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Abstract. We discuss a minimal generalization of the incompressible 
Navier-Stokes equations to describe the solvent flow in an active sus¬ 
pension. To account phenomenologically for the presence of an active 
component driving the ambient fluid flow, we postulate a generic nonlo¬ 
cal extension of the stress-tensor, conceptually similar to those recently 
introduced in granular media flows. Stability and spectral properties of 
the resulting hydrodynamic model are studied both analytically and 
numerically for the two-dimensional (2D) case with periodic boundary 
conditions. Future generalizations of this momentum-conserving theory 
could be useful for quantifying the shear properties of active suspen¬ 
sions. 


1 Introduction 

An active suspension GHa is, roughly speaking, a passive fluid medium that con¬ 
tains at least one ‘micro-swimmer’ species capable of converting chemical into kinetic 
energy. If the swimmer concentration is sufficiently high, their collective dynamics 
can induce rich non-equilibrium flow patterns in the ambient fluid thereby 

causing sig nificant changes in the transport properties and rheological re¬ 

sponse [2ll - l^ of the solvent medium. An intriguing, seemingly generic feature of 
dense active suspensions is the emergence of a characteristic topological defect or 
vortex distance thought to arise from the competition between self¬ 

propulsion, steric and hydrodynamic interactions [^ . Although the microscopic ori¬ 
gins of such dynamical length-scale selection mechanisms are not yet fully understood, 
their experimentally confirmed presence [1-0, lH, Ull suggests that one can effectively 
describe active suspensions in terms of ‘non-local’ higher-than-second-order partial 
differential equations (PDEs) [^, in analogy with well-established continuum theo¬ 
ries of pattern formation in elastic materials , granular media [S^ and convection 
phenomena 

Indeed, recent experimental and theoretical studies ana confirm that a fourth- 
order extension of the Toner-Tu theory [H, Ull can reproduce, both qualitatively 
and quantitatively, many of the main statistical features of dense bacterial sus¬ 
pensions. Specifically, these experiments Bin measured the mean bacterial ve¬ 
locity field u{t,x)^ which can be approximately decomposed in the form u(t,x) ~ 
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v{t, x) + voP{t, x), where v{t, x) is the underlying solvent velocity field and P{t, x) 
denotes the local mean orientation of the bacteria. The parameter vg is the typi¬ 
cal bacterial self-swimming speed relative to the solvent flow (in general, vg is also 
a fluctuatin g q uantity). The bacterial velocity data {it}, obtained by standard PIV 
methods [1,11, 0, were found to agree well with predictions of the incompressible 
fourth-order theory BUI 

V • M = 0 (la) 

{dt + Xgu ■ V)m = —V(p — Aitt^) — /3(m^ — Uq)u + PgV'^u — T 2 (V^)^it, (lb) 

where the pressure p{t, x) is the Lagrange multiplier for the incompressibility (bac¬ 
terial mass conservation) constraint dial) . The parameter Ag describes nematic ad- 
vection and Ai an active pressure contribution. The (/3, uo)-terms correspond to a 
quartic Landau-type velociW potential (3l - l35j| and account for the formation of lo¬ 
cally aligned bacterial jets [ 1 ]. The nonlocal {Fg, r’ 2 )-terms encode passive and active 
stresses due to hydrodynamic and steric interactions, and determine the characteristic 
vortex size Ar — 2Try^r2/{—rg) in the model when 1^2 > 0 and Fg < 0. Conceptually, 
Eqs. m extend the incompressible Toner-Tu theory (33l - l35l| through the additional 
Swift-Hohenberg-type [3l| instability that arises for /g < 0. A similar linear insta¬ 
bility mechanism was recently derived by Grofimann et al. who considered a 
self-propelled particle model with velocity-dependent interaction. With regard to our 
subsequent discussion, it is important to note that Eqs. o, which are defined in the 
rest-frame of the microfluidic channel confining the suspensions, are non-conservative 
due to the aligning /?-term, reflecting the fact that the self-swimming field vgP is not 
a conserved quantity (just as in the classical Toner-Tu model [3^. [3^1. 

Although Eqs. dH) give satisfactory predictions for the bacterial velocity field 
u{t, x) in a stationary setting am, they are of limited use with regard to shear 
experiments, which typically measure the response of the solvent flow in the presence 
of moving boundaries. Aiming to develop a simplified phenomenolo gica l framework 
for the future mathematical description of rheological measurements . we will 

focus here on the complementary problem of constructing effective models for the 
solvent velocity field v(t,x). Specifically, we are interested in identifying a minimal 
extension of the Navier-Stokes (NS) equations that reproduces qualitatively the ex¬ 
perimentally observed, turbulent tracer dynamics in active suspensions [3]. In contrast 
to traditional approaches that build on explicit coup lin gs b etween solvent flow and 
additional orientational order-parameter fields (^lll lla . Hj, [SJ] , we investigate here 
analytically and numerically higher-order ad hoc closure conditions for the stress ten¬ 
sor. Despite some technical differences, our approach shares conceptual similarities 
with the recently proposed, effectively non-local constitutive relations that have led 
to promising progress in the quantitative understanding of porous media flows [S^. 


2 Generalized Navier-Stokes (NS) model 


We focus on a coarse-grained model of active micro-swimmer suspensions, assuming 
that a single velocity field, v{t,x), describes the solvent flow on scales several times 
larger than an individual micro-swimmer. Considering incompressible solvent flow, 
we assume that the dynamics of v[t,x) is governed by the mass and momentum 
conservation laws 


0 = V • 

dtv + {v ■ y)v = —Vp -I- V • (T, 


(2a) 

(2b) 
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with scalar pressure p(t, x) and symmetric stress tensor cr{t, x). As usual, the swim¬ 
mers are assumed to drive solvent flow by modifying the stress field cr. However, 
instead of constructing cr from orientational order-parameter fields [ol-fTll [TsI . [H, , 

we hypothesize that the stress generated by the micro-swimmers can be captured 
through a generic closed-form ansatz 

cT = i:{y,v). ( 3 ) 

In this paper, we will focus on a representative of the class of isotropic traceless 
tensor^ 


i:(V,^) = /(V 2 )[(V^)-f (Vt;)^], ( 4 ) 

where is the Laplace operator, and the scalar function /(•) quantifies the swimmers- 
solvent coupling. Intuitively, Eq. (| 3 ]) can be thought to arise from a truncated gradient 
expansion of an integral kernel representation of the ‘full’ stress tensor cr, similar to a 
Kramers-Moyal expansion In the limit case of a constant function /(V^) = Iq, 
corresponding to a passive isotropic fluid, Eqs. reduce to the standard Navier- 
Stokes equations. 

In the remainder, we will restrict the discussion to symmetric second-order poly¬ 
nomials 


/(v^) = ro-r2(v2) + A(v2)2. (5) 

The constants /q and Ej are assumed to be positive to ensure asymptotic stability, 
whereas the parameter ^2 may have either sign. Nontrivial steady-state flow structures 
emerge for negative values E2 < 0 . Inserting Eqs. ® and ([S]) into Eqs. @, we obtain 
the hydrodynamic equations 


0 = V • V, (6a) 

dtv + (v ■ V)v = -Vp -b - r2V^v + E4VS, (6b) 

where = (V^)" from now on. Below, we analyze the generalized NS equations (O 
on a square domain. Since we aim to understand their bulk behavior, we adopt peri¬ 
odic boundary conditions throughout. 


3 Analytical results 


To obtain some intuition about the generalized NS model (jH]), we first note that its 
linear part supports a stationary vortex lattice of period ^ when E2 is 

negative. This follows from the fact that 


- E2V'‘ -f = 0 


( 7 ) 


\i k = |fc| is one of the roots of Eq -I- E2fc^ -I- = 0 , with real positive roots 

existing only when E2 < 0 . Furthermore, since the nonlinear advective terms on the 
will generally lead to mixing, one can expect to find parameters 
H) produce turbulent mesoscale patterns similar to those observed in 

iS- 


Ihs. of Eq. 
such that Eqs. 
experiments [^, 


^ More generally, one could also consider additional quasi-nematic stress contribu¬ 
tions oc vv — 7|u|'^/d, where I is the d-dimensional identity tensor. Such terms would ef¬ 
fectively rescale the advective derivative and add a kinetic pressure contribution. 
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To perform a more detailed stability analysis, we focus on a periodic square domain 
f2 = [—L/2,L/2]^ and rescale Eqs. dH) by introducing dimensionless quantities 


27r (27r)2ro 

X -)■ —X, t — — — t, 

L 


L 




inTn 


P 


{2n)^r§ 


P, 


k 




One then finds that the dynamics of the model is characterized by the two dimen¬ 
sionless groups 


_ roFi _ ( 271)^02 

and the rescaled Eqs. (O take the form 

V • u = 0 

dtv + [v ■ V)v = —Vp + — 72 V"^u H- 


(9) 


( 10 a) 

( 10 b) 


To ensure stability at short wavelengths, 7 must be always positive. Nontrivial flow 
structures require 72 < 0. We expect that there is a region in the (7, 72 )-parameter 
space that supports a quasi-chaotic steady-state dynamics characterized by the cre¬ 
ation and annihilation of vortices. 

To estimate this parameter region, we investigate how the total kinetic energy, 
E{t) = ^ J^dx |up, varies with time. From the equations of motion with periodic 
boundary conditions, one finds 


E(t) = / dxv ■ dtV 

J a 

—V • -I- pv 

= j da; u • — 72 V"^u-I- 77|V®u) . ( 11 ) 

Jn 

We next insert the Fourier series for the velocity, v{t, x) = k), 

where fc G Z^, to obtain 

E{t) = -(27r)2^A:2(l -f 72^:^ -b77|fc'^)|D(t,fc)p. (12) 

k 

The first and the third term in the brackets are always positive, and therefore dissipate 
energy. If 72 > 0 holds, then E < 0 always; in this case, any solvent flow in the 
system is rapidly damped out. More interestingly, however, when 72 < 0, the active 
component pumps energy into the flow. We may then ask if, at least for some region 
in the (7, 72 )-plane, the energy input and the energy dissipation can balance each 
other. If such a steady-state exists, then the total system energy should fluctuate 
about a constant mean value. More formally, we expect in this case that the mean 
energy change vanishes. 



^ +v ■ (V^u — 72 V^u 4- 772 V®u) 


{E)a,t 


lim lim — 
A—>ooT—^oo Zj 


rT-\-A 


E{t)dt —>• 0, 


(13a) 


and that the time-averaged Fourier coefficients become stationary and isotropic, 

(|*(t,fc)P)AT ^ (|*(fc)P). (13b) 
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Provided that the steady-state is attainable, we have the following energy balance 
equation 

^ fc^(l-I- 72 ^^ + 77 |fc'‘)£(fc) = 0, (14a) 

k 

where the energy spectrum, £(fc), is defined as 

(If dx\v\^\ =^£(A:), (14b) 

\-^Jn / A.T 

yielding 

= ^ E (15) 

h' '.\k'\—k 

Building on the above considerations, we can now analytically estimate the part of 
the ( 7 , 72 )-parameter plane where a steady-state can be reached. Neglecting pressure, 
the linearized version of Eq. (llObI) reads 

dtv = + 77 |V®t;. (16) 

Using the Fourier series for v as before, this is equivalent to 

dtv{t,k) = A{k)v(t,k), (173-) 

where 

A{k) = -fc^(l -I- 72^^ -I- 

=-lllk^ik^-ki){k^-kX). (17b) 

The range of physically reasonable zeros k± can be inferred from stability considera¬ 
tions as follows: 

Fourier modes may become unstable if the micro-swimmers inject a sufficiently 
large amount of energy into the system. In the linearized model, such supercritical 
energy injection leads to divergence at an exponential rate. However, the advective 
term in the full nonlinear system mixes different modes, facilitating energy dissipation 
through the decaying modes. The stability of a given mode is determined by the 
sign of A{k). To observe non-trivial flow structures, we require A{k) > 0 for some 
fc > 0, implying that 72 < 0 and 7 < 1/4, because otherwise (1 -I- 72 / 2 ^ -I- 77 !^^) > 0 
and, hence, no real positive roots k± exist. Furthermore, its is plausible to assume 
that any realistic active system has a long-wavelength cut-off corresponding to the 
largest scale at which energy is collectively injected into the fluid, implying that the 
system should be dissipative in some vicinity of its lowest mode, k = 1. Otherwise, 
there is no room left for an inverse dissipative energy cascade and the energy could 
continuously accumulate at long wavelengths. We therefore demand fc_ > 1, which 
implies that 72 > 3^(1 — i/l — 47 ). Thus, all discrete fc-modes lying in the annular 
region 1 < < /c < inject energy, whereas the complementary set of modes 

dissipates energy. For instance, if we adopt the simplifying assumption that, as in 
classical 2D turbulent flow, £(fc) scales as up to some cut-off mode, which 

we take to be the largest unstable mode, k'^ = (1 + VI ~ 47 ), then, after 

approximating the sum by an integral, Eq. (I14all predicts the critical value of 7 = 0.16. 

In summary, these considerations suggest that an energetically stable, nontrivial 
steady-state can be reached in the parameter range 0.16 < 7 < 0.25 and 0 > 72 > 
3^(1 — 1/1 ~ 47 ). These simple estimates agree well with our numerical results, as 
shown in Fig. [T]4.. 
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4 Stream function formulation and numerical implementation 

To solve Eqs. m numerically, it is convenient to reformulate the dynamical equa¬ 
tions (UM in terms of a stream function. By means of the Helmholtz-Hodge decompo¬ 
sition |4l| . we can express the solvent flow field u as a sum of divergence-free, curl-free 
and harmonic components, 


V = \/<P + \/ A^ + V, (18) 

for some scalar functions and and a harmonic vector field V{t,x) 

satisfying = 0. On periodic domains, harmonic functions are constant, so V 

is interpreted as the fluid center of mass velocity. We will always work in the center 
of mass frame, hence V = 0 from now on. Similarly, <P is also a constant since 
incompressibility implies that = 0, so that ^ is a harmonic function. Thus, 

Eq. (fT51) reduces to v = 'S/ A'F. 

To obtain the evolution equation for 'P, we take the divergence and the curl of 
Eq. (|10b|) . which gives 

(VVlf) ; (VVlf) - = -V^p, (19a) 

dtiV^W) + A (19b) 

The main advantage of this reformulation lies in the fact that the stream function ^ 
is now the only dynamical variable, as the pressure p can always be recovered from ^ 
by solving the Poisson equation (ll9aD . 

The governing equation for the stream function, Eq. CSg, can be solved with 
standard spectral methods [i^. Using the Fourier series representation, the modes 
of W evolve according to 


dtP+ 'N =+^k‘^ ( 20 a) 

where the nonlinear terms are abbreviated by 

l\r = /c-2T{T-i(ifcfc2#) A (20b) 

We integrated Eqs. (uni) with a classical fourth-order Runge-Kutta scheme, and ap¬ 
proximated T(-) by a Discrete Fourier Transform (DFT). The nonlinear term IV is 
evaluated by inverting the DFT, performing the multiplication in position space, and 
then applying the DFT again. The aliasing error generated in this procedure is re¬ 
moved at the expense of using a larger number of Fourier modes and zero-padding 
when necessary (see (4^ for a detailed description of this method; in our simulations, 
we used a symmetric grid of 243 x 243 modes). 


5 Results 

We first investigated the stability of numerical solutions of Eqs. (HOD on a periodic 
L X L square domain by performing systematic scan of the ( 7 , 72 )-parameter plane, 
initializing each simulation run with small random stream function values. In agree¬ 
ment with our analytical considerations in Sec. |3l we found three qualitatively differ¬ 
ent asymptotic behaviors: 

(i) For 72 > 0 or when 7 becomes too large, the solvent dynamics is purely 
dissipative and approaches a stationary state of vanishing flow v = 0 (blue symbols 
in Fig. [TJA). 
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Fig. 1 Simulation results for the rescaled generalized Navier-Stokes model defined in 
Eqs. (1101) . A: Numerical stability analysis in the ( 7 , 72 )-parameter plane. Nontrivial 
energetically stable steady-state solutions exist in the green region (labels D, E and 
F correspond to snapshots shown in panels D-F). The numerical results agree well 
with the analytically predicted region of flow structure formation (red dashed lines) 
obtained in Sec.|31 In the purely dissipative regime (dark blue), the system converges 
to the zero-flow solution. The truncated model becomes unstable (brown region) 
when the dissipation is not able to balance the energy input. B, C: Kinetic energy 
as a function of time, E(<), and corresponding normalized energy spectra, £(fc), for 
the three parameter choices in the bottom row. D-F: Snapshots of the steady-state 
flow stream lines for simulations with periodic boundary conditions. The background 
color represents the associated vorticity fields lj = VAt>, normalized by their maximal 
values. The domain size \s L x L, with L = 600 /rm in physical units. 


(ii) For 72 < 0 and 7 too small, the solutions become unstable, reflected by an 
exponential blow-up of the kinetic energy (brown symbols in Fig. [T]A.). 

(iii) For 72 < 0 and moderate values of 7 the system exhibits quasi-chaotic steady- 
state flow patterns (green symbols in Fig.[T14). The numerically estimated boundaries 
of this physically relevant domain agree well with the analytical estimates (red dashed 
lines) from Sec. [31 

Example snapshots of quantitatively different flow patterns for three parameter 
pairs ( 7 , 72 ) are shown in Fig.[Tp-F. For all three parameter pairs, the kinetic energy 
approaches a constant mean value (Fig.[Tj3). As 7 increases, the corresponding energy 
spectra develop peaks near the linearly most unstable wave number, fc* ^ —( 2772 )“^ 
(black curve F in Fig.jTp). In position space, these peaks correspond to a characteristic 
vortex size ^ k~^ (Fig. |TF). Near 7 = 0.16, the spectra are well approximated by 
Kolmogorov scaling £(fc) oc with sharp cut-off (red curve E in Fig. IDF), in 

agreement with the assumptions made in Sec. [3] to derive the left vertical red dashed 
line in Fig. [T14. 

To relate the dimensionless parameters to physical relevant dimensional values, 
we may fix /q = 10^ /rm^/s and identify the box length with L = 600 /rm. With these 
choices, the typical steady-state speeds are in the range of 1 /nWs to 100 /rm/s, as 
typical of passive tracer particles in dense bacterial suspensions Q- 
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Finally, we would still like to emphasize that the unstable regime (brown symbols 
in Fig. [T]4_) arises from the particular truncated polynomial ansatz in Eq. Future 
quantitative comparison with experiments should focus on reconstructing better ap¬ 
proximations of the function / or, more generally, cr. Conversely, however, stability 
criteria provide useful physical constraints for effective models that can be utilized in 
parameter estimation procedures. 


6 Summary 

We proposed and analyzed a minimal generalization of the Navier-Stokes equation to 
describe the solvent flow in an active suspension. The main assumption underlying this 
simple momentum-conserving model is that, on scales larger than the swimming cells 
or filaments, the complex fluid-swimmer interactions can be effectively captured by a 
generalized form of the stress energy tensor, which can be expanded in terms of higher- 
order differential operators. In this contribution, we focussed on a simple example, 
corresponding to a sixth-order PDE, that produces turbulent flow features that are 
qualitatively similar to those observed through passive tracer-particle tracking in 
recent experiments Q- 

A future goal is to embed this class of models into a shear-flow setting as rele¬ 
vant to viscosity measurements pll . . This problem is conceptually and numerically 
nontrivial as appropriate boundary conditions need to be identified and implemented. 
Notwithstanding, such a phenomenological approach may help us to progress towards 
a sufficiently-general-yet-reasonably-simple framework for the classification of rheo¬ 
logical observations in active suspensions. Erom a practical perspective, an interesting 
challenge will be to reconstruct an empirical form of the generalized stress-tensor cr 
from experimentally measured solvent flow data. 
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